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Structure of cold dense matter at subnuclear densities is investigated by quantum molecular 
dynamics (QMD) simulations. We succeeded in showing that the phases with slab-like and rod-like 
nuclei etc. can be formed dynamically from hot uniform nuclear matter without any assumptions 
on nuclear shape. We also observe intermediate phases, which has complicated nuclear shapes. 
Geometrical structures of matter are analyzed with Minkowski functionals, and it is found out that 
intermediate phases can be characterized as ones with negative Euler characteristic. Our result 
suggests the existence of these kinds of phases in addition to the simple "pasta" phases in neutron 
star crusts. 
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Collapse-driven supernovae (SNe) and the following 
formation of neutron stars (NSs) are the most dramatic 
processes during the stellar evolution. These objects pro- 
vide not only astrophysically significant phenomena but 
also interesting material phases inside them; both are 
strongly connected with each other. Investigating the 
properties of nuclear matter under extreme condition is 
one of the essential subjects to clarify the mechanism of 
collapse-driven SNe [Q and the structure of NS crusts . 
This subject is also interesting in terms of the fundamen- 
tal problem of the complex fluids of nucleons. 

At subnuclear densities, where nuclei are about to melt 
into uniform matter, it is expected that the energeti- 
cally favorable configuration possesses interesting spatial 
structures such as rod-like and slab-like nuclei and rod- 
like and spherical bubbles etc., which are called nuclear 
"pasta" . This picture was first proposed by Ravenhall 
et al. H and Hashimoto et al. These works, which 
are based on free energy calculations with liquid drop 
models assuming some specific nuclear shapes, clarify 
that the most energetically stable nuclear shape is deter- 
mined by a subtle balance between nuclear surface and 
Coulomb energies. Although detailed aspects of phase 
diagrams vary with nuclear models (s), the realization 
of the "pasta" phases as energy minimum states can be 
seen in a wide range of nuclear models and basic feature 
of phase diagrams is universal || [h]], i.e. with increasing 
density, the shape of the nuclear matter region changes 
like sphere — > cylinder — > slab — > cylindrical hole — > 
spherical hole — > uniform. This feature is also reproduced 
by Thomas-Fermi calculations [|[ 0, [§• 

The phases with these exotic nuclear structures, if they 
are realized in NS crusts or SN cores, bring about many 
astrophysical consequences. As for those to NS phenom- 
ena, it is interesting to note the relevance of nonspheri- 
cal nuclei in neutron star matter (NSM) to pulsar glitches 
[p| and cooling of NSs ||. For supernova matter (SNM), 



"pasta" phases are expected to affect the neutrino trans- 
port H |l(| and hydrodynamics in SN cores J5J. 

Though the properties of "pasta" phases in equilibrium 
state have been investigated actively, the formation and 
the melting processes of them have not been discussed 
except for some limited cases based on perturbative ap- 
proaches |^|, O]. It is important to adopt a microscopic 
and dynamical approach which allows arbitrary nuclear 
structures to understand these processes of nonsphcri- 
cal nuclei. At finite temperatures, it is considered that 
not only nuclear surface becomes obscure but also nu- 
clei of various shapes may coexist. Therefore, it is nec- 
essary to incorporate density fluctuations without any 
assumptions on nuclear shape to investigate the proper- 
ties of "pasta" phases at finite temperatures. Although 
some previous works Q do not assume nuclear struc- 
ture, they can not incorporate fluctuations of nucleon 
distributions satisfactorily because they are based on the 
Thomas- Fermi calculation which is one-body approxima- 
tion. In addition, only one structure is contained in the 
simulation box in these works, there are thus possibilities 
that nuclear shape is strongly affected by boundary effect 
and some structures are prohibited implicitly. 

We have started studying dense matter at subnuclear 
densities by quantum molecular dynamics (QMD) (ilf| 
which is one of the molecular dynamics (MD) approaches 
for nucleon many body systems (see, e.g., Ref. for 
review). MD for nucleons including QMD, which is a 
microscopic and dynamical method without any assump- 
tions on nuclear structure, is suitable for incorporating 
fluctuations of particle distributions. The final aim of 
our study is to understand the above mentioned forma- 
tion and melting processes of nonspherical nuclei and to 
investigate the properties of matter consists of nonspher- 
ical nuclei at finite temperatures. For the first step, the 
question posed in this paper is as follows: Can the phases 
with nonspherical nuclei be formed dynamically? 
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There is a pioneering work by Maruyama et al. |Q 
which attempts to investigate the structure of matter at 
subnuclear densities by QMD. Unfortunately, they did 
not treat the Coulomb interaction consistently and could 
not anneal the system successfully. As a result, they 
could not reproduce the phases with nonspherical nuclei 
of simple structures. In the present work, we improve 
the above shortcomings and obtain phase diagrams for 
cold dense matter of the proton fraction x = 0.3 and 0.5 
at subnuclear densities. We also discuss nuclear shape 
changes with morphological measures. 

While there are a lot of versions of MD for fermions, 
we choose QMD among them based on trade-off between 
calculational amounts and accuracies. The typical length 
scale I of inter-structure is I ~ 10 fm and the density 
region of interest is just below the normal nuclear density 
po = 0.165 fm -3 . The required nucleon number N in 
order to reproduce n structures in the simulation box is 
about N ~ po(nl) 3 (for slabs), it is thus desirable that 
we prepare nucleons of order 10000 if we try to reduce 
boundary effects. While it is very hard task to treat such 
a large system by, for example, FMD and AMD (see, 
e.g., Ref. G3] and references therein) whose calculational 
amounts scale as ~ iV 4 , it is feasible to do it by QMD 
whose calculational amounts scale as ~ iV 2 . It is also 
noted that we mainly focus on macroscopic structures; 
the exchange effect would not be so important for them. 
Therefore, QMD, which is less elaborate in treating the 
exchange effect, has the advantage of other models. 

We have performed QMD simulations of infinite (n,p,e) 
system with fixed proton fraction x = 0.3 and 0.5 for var- 
ious nucleon densities p ( the density region is 0.05 - 1.0 
po)- We set a cubic box which is imposed periodic bound- 
ary conditions in which 2048 nucleons (1372 nucleons in 
some cases) are contained. The relativistic degenerate 
electrons which ensure the charge neutrality are treated 
as a uniform background and the Coulomb interaction 
is calculated by the Ewald method (see, e.g., Ref. 16 ) 
which enables us to sum contributions of long-range in- 
teractions in a system with periodic boundary conditions. 
For nuclear interaction, we use an effective Hamiltonian 
developed by Maruyama et al. (medium EOS model) 
fl2"f which reproduce the bulk properties of nuclear mat- 
ter and the properties of stable nuclei, especially heavier 
ones, i.e. binding energy and root-mean-square radius. 

We first prepare an uniform hot nucleon gas at ksT ~ 
20 MeV as an initial condition which is equilibrated for 
~ 500 — 2000 fm/c in advance. In order to realize the 
ground state of matter, wc then cool it down slowly for 
O(10 3 — 10 4 ) fm/c keeping the nucleon density constant 
by frictional relaxation method etc. until the temper- 
ature gets ~ 0.1 MeV or less. Note that any artificial 
fluctuations are not given during the simulation. 

The QMD equations of motion with friction terms are 
solved using the fourth-order Gear predictor-corrector 
method in conjunction with multiple time step algorithm 



fl6|| . Integration time steps At are set to be adaptive in 
the range of At < 0.1 — 0.2 fm/c depending on the degree 
of convergence. At each step, the correcting operation is 
iterated until the error of position Ar and the relative 
error of momentum Ap/p get smaller than 10 -6 , where 
Ar and Ap/p are estimated as the maximum values of 
correction among all particles. Computer systems which 
we use are equipped with MD-GRAPE II. 

Shown in Figs. 1 and 2 are the resultant nucleon dis- 
tributions of cold matter at x = 0.5 and 0.3, respec- 
tively. We can see from these figures that the phases 
with rod-like and slab-like nuclei, cylindrical and spheri- 
cal bubbles, in addition to the one with spherical nuclei 
are reproduced in the both cases. We here would like 
to mention some reasons of discrepancies between the 
present result and the result obtained by Maruyama et al. 
which says "the nuclear shape may not have these simple 
symmetries" jl2). The most crucial reason seems to be 
the difference in treatment of the Coulomb interaction. 
In the present simulation, we calculate the long range 
Coulomb interaction in a consistent way using the Ewald 
method. For the system of interest where the Thomas- 
Fermi screening length is comparable to or larger than 
the size of nuclei, this treatment is more adequate than 
that with introducing an artificial cutoff distance as in 
Ref. |Q. The second reason would be the difference 
in the relaxation time scales r. In our simulation, we 
can reproduce the bubble-phases (see d and e of Figs. 1 
and 2) with r ~ 10 3 fm/c and the nucleus-phases (see 
b and c of Figs. 1 and 2) with r ~ O(10 4 ) fm/c. How- 
ever, the matter in the density region corresponding to a 
nucleus-phase is quenched in an amorphous state when 
r < 10 3 fm/c. In the present work, we take r much larger 
than typical time scale r t h ~ O(100) fm/c for nucleons 
to thermally diffuse the distance of I ~ 10 fm at p ~ po 
and ksT ~ 1 MeV. This temperature is below the typ- 
ical value of the liquid-gas phase transition temperature 
in the density region of interest, it is thus considered that 
our results are thermally relaxed in a satisfying level. 

Phase diagrams of matter in the ground state are 
shown in Figs. 3 (a) and (b) for x = 0.5 and 0.3, respec- 
tively. As can be seen from these figures, the obtained 
phase diagrams basically reproduce the sequence of the 
energetically favored nuclear shapes predicted by simple 
discussions j| which only take account of the Coulomb 
and surface effects; this prediction is that the nuclear 
shape changes like sphere — > cylinder — > slab — > cylin- 
drical hole — > spherical hole — > uniform, with increasing 
density. Comparing Figs. 3 (a) and (b), we can see that 
the phase diagram shifts towards the lower density side 
with decreasing x, which is due to the tendency that as 
the nuclear matter at larger neutron excess, the satura- 
tion density is lowered. It is remarkable that the density 
dependence of the nuclear shape, except for cylindrical 
bubbles (just in the case of x = 0.3) and spherical nuclei 
and bubbles, is quite sensitive and phases with interme- 
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FIG. 1: The nucleon distributions of typical phases with simple structures of cold matter at x — 0.5; (a) sphere phase, O.lpo 
(D = 43.65 fm, TV = 1372); (b) cylinder phase, 0.18p (D = 41.01 fm, TV = 2048); (c) slab phase, 0.4p (D = 31.42 fm, 
TV = 2048); (d) cylindrical hole phase, 0.5p (D = 29.17 fm, TV = 2048) and (e) spherical hole phase, 0.6p (D = 27.45 fm, 
TV = 2048), where D is the box size. The red particles represent protons and the green ones represent neutrons. 




FIG. 2: Same as Fig. 1 at x = 0.3; (a) sphere phase, O.lpo {D = 49.88 fm, TV = 2048); (b) cylinder phase, 0.18p (D = 43.58 fm, 
TV = 2048); (c) slab phase, 0.35p (D = 32.85 fm, TV = 2048); (d) cylindrical hole phase, 0.5p (D = 29.17 fm, TV = 2048) and 
(e) spherical hole phase, 0.55po (D = 28.26 fm, TV = 2048). The red particles represent protons and the green ones represent 
neutrons. 
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FIG. 3: Phase diagrams of cold matter at x — 0.5 (a) and x = 
0.3 (b). Matter is unstable against phase separation in the 
density region shown as kt < 0, where kt is the isothermal 
compressibility. The symbols SP, C, S, CH and SH stand 
for nuclear shapes, i.e. sphere, cylinder, slab, cylindrical hole 
and spherical hole, respectively. The parentheses (A,B) show 
intermediate phases between A-phase and B-phase suggested 
in this work. They have complicated structures different from 
those of both A-phase and B-phase. Simulations have been 
carried out at densities denoted by small circles. 



diate nuclear shapes which are not simple as shown in 
Figs. 1 and 2 are observed in two density regions: one is 
between the cylinder phase and the slab phase, the other 
is between the slab phase and the cylindrical hole phase. 
We note that these phases are different from coexistence 



phases with nuclei of simple shapes and we will referred 
to them as "intermediate phases" . 

To extract the morphological characteristics of the nu- 
clear shape changes and the intermediate phases, we in- 
troduce the Minkowski functionals (see, e.g., Ref. |17j 
and references therein) as geometrical and topological 
measures of the nuclear surface. Let us consider a ho- 
mogeneous body K G 7Z in the d-dimensional Eucledian 
space, where 1Z is the class of such bodies. Morphologi- 
cal measures are defined as functionals tp : 1Z — ► R which 
satisfy the following three properties: (1) Motion invari- 
ance, i.e. tp(K) — tp(gK), where g denotes any trans- 
lations and rotations. (2) Additivity, i.e. tp(Ki U K-i) = 
ip(Ki)+ip(K 2 )-ip{KinK 2 ), where Ki,K 2 G K. (3) Con- 
tinuity, i.e. lim n -*.ao<p{Kn) = v(^) if limn^^Kn = K, 
where if is a convex body and {K n } is a sequence of 
convex bodies. Hadwiger's theorem states that there are 
just d+1 independent functionals which satisfy the above 
properties; they are known as Minkowski functionals. In 
three dimensional space, four Minkowski functionals are 
related to the volume, the surface area, the integral mean 
curvature and the Euler characteristic. 

Here, we particularly focus on the integral mean cur- 
vature and the Euler characteristic; the results of other 
quantities will be discussed elsewhere. Both are de- 
scribed by surface integrals of the following local quan- 
tities, the mean curvature H — («i + K2)/2 and the 
Gaussian curvature G = k±K2, i.e. J gK HdA and x = 
5tF Ibk GdA, where Ki and k 2 are the principal curva- 
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tures and dA is the area element of the surface of the 
body K. The Euler characteristic \ is a purely topolog- 
ical quantity and 

X = (number of isolated regions) — (number of tunnels) 
+ (number of cavities). (1) 

Thus x > for the sphere and the spherical hole phases 
and the coexistence phase of spheres and cylinders, and 
X = for the other ideal "pasta" phases, i.e. the cylin- 
der, the slab and the cylindrical hole phases. We intro- 
duce the area-averaged mean curvature, (H) = \j HdA, 
and the Euler characteristic density, x/V, as normalized 
quantities, where V is the volume of the whole space. 

We calculate these quantities by the following proce- 
dure. We first construct proton and nucleon density dis- 

2 2 

tributions p p (r) — | <3> p (r) | and p(r) = |3>(r)| from data 
of the centers of position of the nucleons, where $ p (r) and 
$(r) are the QMD trial wave functions of protons and nu- 
cleons (see Ref. Q). We set a threshold proton density 
p Pith and then calculate /(p P ,th) = V r (p Pjt h)M(pp,th), 
where V(p Pl th) and j4(p p ,th) are the volume and the sur- 
face area of the regions in which p p (r) > p p . t h- We 

find out the value p Pjt h = p£, t h where jpr^f(Pp,th) = 
and define the regions in which p p (r) > p* th as nuclear 
regions. For spherical nuclei, for example, p* th corre- 
sponds to a point of inflection of a radial density distri- 
bution. In the almost whole phase-separating region, the 
values of p* h distribute in the range of about 0.7 — 0.9po 
in the both cases of a; = 0.5 and 0.3, where p* h is the 
threshold nucleon density corresponds to p* th . We then 
calculate A, J HdA and x f° r the determined nuclear 
surface. We evaluate A by the triangle decomposition 
method, J HdA by the algorithm shown in Ref. in 
conjunction with a calibration by correction of surface 
area, and x by the algorithm of Ref. jlTj and by that 
of counting deficit angles |D| which are confirmed that 
both of them give the same results. 

We have plotted the obtained p dependence of (H) 
and of x/V f° r the surface of p t h = p* th in Fig. 4. In 
addition to the values of (H) for the surface of p t h = Pth' 
we have also investigated those for the surface of pth = 
p* th ± 0.05po to examine the extent of the uncertainties 
of this quantity which stem from the arbitrariness in the 
definition of the nuclear surface, but we could not observe 
remarkable differences from the values for p t h = Pth (they 
were smaller than 0.015/fm). We could not see these 
kinds of uncertainties in x/V except for the densities near 
below the density at which matter turns into uniform. 

The behavior of (H) shows that it decreases almost 
monotonously from positive to negative with increasing 
p until the matter turns into uniform. The densities cor- 
respond to (H) ~ are about 0.4 and 0.35po for x = 0.5 
and 0.3, respectively; these values are consistent with 
the density regions of the phase with slab-like nuclei (see 
Fig. 3). As mentioned previously, x/V is actually pos- 
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FIG. 4: Density dependence of the Minkowski functionals of 
cold matter at x = 0.5 (a and c) and x = 0.3 (b and d). 



itive in the density regions corresponding to the phases 
with spherical nuclei, coexistence of spherical and cylin- 
drical nuclei, and spherical holes because of the existence 
of isolated regions. As for those corresponding to the 
phases with cylindrical nuclei, planar nuclei and cylin- 
drical holes, x/V — 0. The fact that the values of x/V 
are not exactly zero for nucleon distributions shown as 
slab phases in Figs. 1 and 2 reflects the imperfection 
of these "slabs" , which is due to the small nuclear parts 
which connect the neighboring slabs. However, we can 
say that the behavior of x/V depicted in Figs. 4(c) and 
4(d) shows that x/V is negative in the density regions of 
the intermediate phases, even if we take account of the 
imperfection of the obtained nuclear shapes and the un- 
certainties of the definition of the nuclear surface. This 
means that the intermediate phases consist of nuclear 
surfaces which are saddle-like at each point on average 
and they consist of each highly connected nuclear and 
gas regions due to a lot of tunnels (see Eq. ([j])). 

Let us now refer to discrepancies from the results of 
previous works which do not assume nuclear structure 
H the intermediate phases can not be seen in these 
works. We can give following two reasons for the discrep- 
ancies: (1) These previous calculations are based on the 
Thomas-Fermi approximation which can not sufficiently 
incorporate fluctuations of nucleon distributions. This 
shortcoming may result in favoring nuclei of smoothed 
simple shapes than in the real situation. (2) There is a 
large possibility that some highly connected structures 
which have two or more substructures in a period are 
neglected in these works because only one structure is 
contained in a simulation box. 

If the phases with highly connected nuclear and bub- 
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ble regions are realized as the most energetically stable 
state, we can say that it is not unnatural thing JTs(] . It 
is considered that, for example, a phase with perforated 
slab-like nuclei, which has negative x/V , could be more 
energetically stable than that with extremely thin slab- 
like nuclei. The thin planar nucleus costs surface-surface 
energy which stems from the fact that nucleons in it feel 
surfaces of both sides. We have to examine the existence 
of the intermediate phases by more extensive simulations 
with larger nucleon numbers and with longer relaxation 
time scales in the future. 

Here we would like to discuss astrophysical conse- 
quences of our results. Pethick and Potekhin have 
pointed out that "pasta" phases with rod-like and slab- 
like nuclei are analogous to the liquid crystals according 
to the similarity of the geometrical structures (ll| . It can 
also be said that the intermediate phases observed in the 
present work are "sponge" -like phases because they have 
both highly connected nuclear and bubble regions shown 
as x/V < 0. The elastic properties of the sponge-like in- 
termediate phases are qualitatively different from those 
of the liquid crystal-like "pasta" phases because the for- 
mer ones do not have any directions in which restoring 
force does not act, on the other hand the latter ones have. 
Our results suggest that the intermediate phases occupy 
a significant fraction of the density region in which non- 
spherical nuclei can be seen (see Fig. 3). If this is also 
true for more neutron- rich case as x ~ 0.1, it leads to 
increasing of the maximum clastic energy that can be 
stored in the NS crust than that in the case of the all 
nonspherical nuclei have simple structures. Besides, the 
cylinder and the slab phases which are liquid crystal-like 
lie between the sponge-like intermediate phases or the 
crystalline solid-like phase, and the releasing of the strain 
energy would, in consequence, concentrate in the domain 
of these liquid crystal-like phases. The above mentioned 
effects of the intermediate phases should be taken into 
account in considering the crust dynamics of starquakes 
etc. if these phases exist in NSM. In the context of pulsar 
glitch phenomena, the effects of the sponge-like nuclei on 
pinning rate of superfluid neutron vortices also have yet 
to be investigated. 

For neutrino cooling of NSs, some version of the direct 
URCA process which is suggested by Lorenz et al. || 
that this might be allowed in the "pasta" phases would 
be suppressed in the intermediate phases. This is due to 
the fact that the proton spectrum is no longer continu- 
ous in the sponge-like nuclei. The last point which we 
would like to mention is about the effects of the inter- 
mediate phases on neutrino trapping in SN cores. The 
nuclear parts connect over a wide region which is much 
larger than that characterized by the typical neutrino 
wave length ~ 20 fm. Thus the neutrino scattering pro- 
cesses are no longer coherent in contrast to the case of the 
spherical nuclei, and this may, in consequence, reduce the 
diffusion time scale of neutrinos as in the case of "pasta" 



phases with simple structures. This reduction softens the 
SNM and would thus act to enhance the amount of the 
released gravitational energy. 

Our calculations demonstrate that the "pasta" phases 
can be formed dynamically from hot uniform matter in 
the proton-rich cases of x = 0.5 and 0.3 without any 
assumptions on nuclear shape. This suggests that the 
existence of these phases in NS crusts because they cool 
down keeping the local thermal equilibrium after proto- 
NSs are formed and their cooling time scale is much larger 
than the relaxation time scale of our simulations. This 
conclusion has to be confirmed in more neutron-rich cases 
of x ~ 0.1 in the future. Our results also suggest that 
the existence of the highly connected intermediate phases 
which are characterized as x/V < 0- This provides a 
vivid picture that NS inner crusts which consist of dense 
matter at subnuclear densities may be rich in properties 
due to the possibilities of a variety of material phases. 
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